/****
   Copyright 2005-2007, Moshe Looks and Novamente LLC

   Licensed under the Apache License, Version 2.0 (the "License");
   you may not use this file except in compliance with the License.
   You may obtain a copy of the License at

       http://www.apache.org/licenses/LICENSE-2.0

   Unless required by applicable law or agreed to in writing, software
   distributed under the License is distributed on an "AS IS" BASIS,
   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   See the License for the specific language governing permissions and
   limitations under the License.
****/

#include "reduct/mixed_rules.h"
#include "combo/tree_type.h"
#include <map>

namespace reduct {
  typedef vtree::sibling_iterator sib_it;
  typedef vtree::iterator pre_it;

  //0<c*x -> 0<x if c>0
  //0<c*x -> 0<-1*x if c<0
  //WARNING : this rule is deprecated, use reduce_gt_zero_prod instead
  void reduce_gt_zero_times_const::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::times) {
	for(sib_it sib = it_child.begin(); sib != it_child.end();) {
	  if(is_contin(*sib)) {
	    contin_t c = get_contin(*sib);
	    if(c>0) {
	      if(it_child.has_one_child()) {
		tr.erase_children(it);
		*it=id::logical_true;
		return;
	      }
	      else sib = tr.erase(sib);
	    }
	    else if(c<0) {
	      if(it_child.has_one_child()) {
		tr.erase_children(it);
		*it=id::logical_false;
		return;
	      }
	      else {
		*sib=-1.0;
		++sib;
	      }
	    }
	    else {//that is c==0
	      tr.erase_children(it);
	      *it=id::logical_false;
	      return;
	    }
	  }
	  else ++sib;
	}
      }
    }
  }

  //0<c*x*x -> false if c<0
  //or more generally 0<c*x^p_x*y^p_y*exp(... where p_x, p_y... are divisible by 2
  void reduce_gt_zero_pair_power::operator()(vtree& tr,vtree::iterator it) const {
    //associate for each subtree its parity
    typedef std::map<pre_it, bool,
      util::lexicographic_subtree_order<vertex> > subtree_parity;
    typedef subtree_parity::iterator subtree_parity_it;
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      bool is_neg_const = false;
      subtree_parity sp;
      if(*it_child==id::times) {
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib) {
	  if(is_contin(*sib)) {
	    if(get_contin(*sib)<0)
	      is_neg_const = !is_neg_const;
	  }
	  else {
	    pre_it pre_sib = sib;
	    subtree_parity_it sp_it = sp.find(pre_sib);
	    if(sp_it==sp.end()) {
	      std::pair<pre_it, bool> p(pre_sib, false);
	      sp.insert(p);
	    }
	    else {
	      bool& par = sp_it->second;
	      par = !par;
	    }
	  }
	}
	if(is_neg_const) {
	  if(!sp.empty()) {
	    for(subtree_parity_it sp_it=sp.begin();sp_it!=sp.end();++sp_it) {
	      if(!sp_it->second) return;
	    }
	  }
	  //if still in the code then all sp_it are positive
	  //and so the expression is negative
	  tr.erase_children(it);
	  *it=id::logical_false;
	}
      }
    }
  }

  //0<c/x -> 0<x if c>0
  //0<c/x -> 0<-1*x if c<0
  //WARNING : this rule is deprecated, use reduce_gt_zero_div instead
  void reduce_gt_zero_const_div::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::div) {
	assert(it_child.number_of_children()==2);
	pre_it num = it_child.begin();
	if(is_contin(*num)) {
	  contin_t c = get_contin(*num);
	  if(c>0) {
	    pre_it denom = tr.child(it_child, 1);
	    tr.move_before(it_child, denom);
	    tr.erase(it_child);
	  }
	  else if(c<0) {
	    *num = -1.0;
	    *it_child = id::times;
	  }
	  else {//that is c==0
	    tr.erase_children(it);
	    *it=id::logical_false;
	    return;
	  }
	}
      }
    }
  }

  //0<log(x) -> 0<-1+x
  void reduce_gt_zero_log::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::log) {
	assert(it_child.has_one_child());
	*it_child = id::plus;
	tr.append_child(it_child, -1.0);
      }
    }
  }

  //0<exp(x) -> true
  void reduce_gt_zero_exp::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::exp) {
	assert(it_child.has_one_child());
	tr.erase_children(it);
	*it = id::logical_true;
      }
    }
  }

  //0<-1*exp(x) -> false
  //WARNING : this rule is deprecated, use reduce_gt_zero_prod instead
  void reduce_gt_zero_minus_exp::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::times) {
	bool is_neg_const = false;
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib) {
	  if(is_contin(*sib)) {
	    if(get_contin(*sib) < 0)
	      is_neg_const = !is_neg_const;
	  }
	  else if(*sib!=id::exp)
	    return;
	}
	if(is_neg_const) {
	  tr.erase_children(it);
	  *it = id::logical_false;
	}
      }
    }
  }

  //0<y*exp(x) -> 0<y
  //WARNING : this rule is deprecated, use reduce_gt_zero_prod instead
  void reduce_gt_zero_prod_exp::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::times) {
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib) {
	  if(*sib==id::exp)
	    tr.erase(sib);
	}
	if(it_child.is_childless()) {
	  tr.erase_children(it);
	  *it = id::logical_true;
	}
      }
    }
  }

  //0<c+sin(x) -> true if c>1
  //0<c+sin(x) -> false if c<=1
  //WARNING : this rule is deprecated, use reduce_gt_zero_sum_sin instead
  void reduce_gt_zero_const_sum_sin::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::plus) {
	contin_t total = 0.0;
	int sin_count = 0;
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib) {
	  if(is_contin(*sib))
	    total = total + get_contin(*sib);
	  else if(*sib==id::sin)
	    ++sin_count;
	  else return;
	}
	if(sin_count > 0) {
	  if(total > (contin_t)sin_count) {
	    tr.erase_children(it);
	    *it = id::logical_true;
	  }
	  else if(total <= (contin_t)-sin_count) {
	    tr.erase_children(it);
	    *it = id::logical_false;
	  }
	}
      }
    }
  }

  //0<impulse(x) -> x
  void reduce_gt_zero_impulse::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::impulse) {
	assert(it_child.has_one_child());
	pre_it new_it = it_child.begin();
	*it = *new_it;
	tr.reparent(it, new_it);
	tr.erase(it_child);
      }
    }
  }

  //prod impulse(x_i)^p_i -> prod impulse(x_i)
  void reduce_impulse_power::operator()(vtree& tr,vtree::iterator it) const {
    typedef std::set<pre_it, util::lexicographic_subtree_order<vertex>
      > subtree_quotient;
    if(*it==id::times) {
      subtree_quotient sq;
      for(sib_it sib = it.begin(); sib != it.end();) {
	if(*sib==id::impulse) {
	  assert(sib.has_one_child());
	  if(sq.count(pre_it(sib)) > 0) {
	    sib = tr.erase(sib);
	  }
	  else {
	    sq.insert(pre_it(sib));
	    ++sib;
	  }
	}
	else ++sib;
      }
      if(it.has_one_child()) {
	*it = *it.begin();
	tr.erase(tr.flatten(it.begin()));
      }
    }
  }

  //prod impulse(x_i) * z -> impulse(and x_i) * z
  void reduce_impulse_prod::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::times) {
      pre_it fic; //child of the first impulse
      pre_it new_and;
      for(sib_it sib = it.begin(); sib != it.end();) {
	if(*sib==id::impulse) {
	  assert(sib.has_one_child());
	  if(tr.is_valid(fic)) {
	    if(!tr.is_valid(new_and))
	      new_and = tr.wrap(fic, id::logical_and);
	    tr.move_after(fic, pre_it(sib.begin()));
	    sib = tr.erase(sib);
	  }
	  else {
	    fic = pre_it(sib.begin());
	    ++sib;
	  }
	}
	else ++sib;
      }
      if(it.has_one_child()) {
	*it = *it.begin();
	tr.erase(tr.flatten(it.begin()));
      }
    }
  }

  //0<(sum impulse(x_i)) -> or x_i
  void reduce_impulse_sum::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::plus) {
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib)
	  if(*sib!=id::impulse) return;
	//we are now sure that all arguments of + are impulses
	*it = id::logical_or;
	tr.erase(tr.flatten(it_child));
	for(sib_it sib = it.begin(); sib != it.end(); ++sib) {
	  sib = tr.erase(tr.flatten(sib));
	}
      }
    }
  }

  //if(x 1 0) -> impulse(x)
  //if(x 0 1) -> impulse(NOT(x))
  void reduce_contin_if_to_impulse::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it cond = tr.child(it, 0);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      if(is_contin(*b1) && is_contin(*b2)) {
	contin_t cb1 = get_contin(*b1);
	contin_t cb2 = get_contin(*b2);
	if(cb1==1.0 && cb2==0.0) {
	  *it = id::impulse;
	  tr.erase(b1);
	  tr.erase(b2);
	}
	else if(cb1==0.0 && cb2==1.0) {
	  tr.wrap(cond, id::logical_not);
	  *it = id::impulse;
	  tr.erase(b1);
	  tr.erase(b2);
	}
      }
    }
  }

  //if(true x y) -> x
  //if(false x y) -> y
  //if(x if(x y z) w) -> if(x y w)
  //if(x y if(x z w)) -> if(x y w)
  void reduce_contin_if::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it cond = tr.child(it, 0);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      if(*cond==id::logical_true) {
	*it = *b1;
	tr.flatten(b1);
	tr.erase(cond);
	tr.erase(b1);
	tr.erase(b2);
      }
      else if(*cond==id::logical_false) {
	*it = *b2;
	tr.flatten(b2);
	tr.erase(cond);
	tr.erase(b1);
	tr.erase(b2);
      }
      else if(*b1==id::contin_if) {
	assert(tr.number_of_children(b1)==3);
	pre_it b1_cond = tr.child(b1, 0);
	pre_it b1_b1 = tr.child(b1, 1);
	if(tr.equal_subtree(cond, b1_cond)) {
	  tr.move_before(b1, b1_b1);
	  tr.erase(b1);
	}
      }
      else if(*b2==id::contin_if) {
	assert(tr.number_of_children(b2)==3);
	pre_it b2_cond = tr.child(b2, 0);
	pre_it b2_b2 = tr.child(b2, 2);
	if(tr.equal_subtree(cond, b2_cond)) {
	  tr.move_before(b2, b2_b2);
	  tr.erase(b2);
	}
      }
    }
  }
  
  //op(if(x y1 z1) if(x y2 z2)) -> if(x op(y1 y2) op(z1 z2))
  void reduce_op_contin_if::operator()(vtree& tr,vtree::iterator it) const {
    typedef std::set<pre_it, util::lexicographic_subtree_order<vertex>
      > subtree_quotient;
    typedef subtree_quotient::iterator subtree_quotient_it;
    if(*it==id::div) {
      assert(tr.number_of_children(it)==2);
      pre_it num = tr.child(it, 0);
      pre_it denum = tr.child(it, 1);
      if(*num==id::contin_if && *denum==id::contin_if) {
	assert(tr.number_of_children(num)==3);
	assert(tr.number_of_children(denum)==3);
	pre_it num_cond = tr.child(num, 0);
	pre_it denum_cond = tr.child(denum, 0);
	if(tr.equal_subtree(num_cond, denum_cond)) {
	  *it=id::contin_if;
	  //move x as first child of the new if
	  tr.move_before(num, num_cond);
	  tr.erase(denum_cond);
	  *num=id::div;
	  *denum=id::div;
	  tr.swap(pre_it(tr.child(num, 1)), pre_it(tr.child(denum, 0)));
	}
      }
    }
    else if(*it==id::plus || *it==id::times) {
      subtree_quotient sq;
      for(sib_it sib = it.begin(); sib != it.end(); ) {
	if(*sib==id::contin_if) {
	  assert(tr.number_of_children(sib)==3);
	  pre_it sib_cond = tr.child(sib, 0);
	  subtree_quotient_it sqi = sq.find(sib_cond);
	  if(sqi==sq.end()) {
	    sq.insert(pre_it(sib_cond));
	    ++sib;
	  }
	  else {
	    pre_it pivot_cond = *sqi;
	    pre_it pivot_b1 = tr.next_sibling(pivot_cond);
	    pre_it pivot_b2 = tr.next_sibling(pivot_b1);
	    tr.erase(tr.child(sib, 0));
	    builtin it_builtin = get_builtin(*it);
	    if(*pivot_b1==it_builtin)
	      pivot_b1 = pivot_b1.begin();
	    else tr.wrap(pivot_b1, *it);
	    tr.move_before(pivot_b1, pre_it(tr.child(sib, 0)));
	    if(*pivot_b2==it_builtin)
	      pivot_b2 = pivot_b2.begin();
	    else tr.wrap(pivot_b2, *it);
	    tr.move_before(pivot_b2, pre_it(tr.child(sib, 0)));
	    sib = tr.erase(sib);
	  }
	}
	else ++sib;
      }
      if(it.has_one_child()) {
	*it=*it.begin();
	tr.erase(tr.flatten(it.begin()));
      }
    }
  }

  //contin_if(x y y) -> y
  void reduce_contin_if_equal_branch::operator()(vtree& tr,
						 vtree::iterator it) const {
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it cond = tr.child(it, 0);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      if(tr.equal_subtree(b1, b2)) {
	*it=*b1;
	tr.erase(tr.flatten(b1));
	tr.erase(cond);
	tr.erase(b2);
      }
    }
  }

  //contin_if(x op(y z) op(y w)) -> op(y contin_if(x z w))
  //op in {+, *, /}. If op is / the order of argument is respected
  void reduce_contin_if_inner_op::operator()(vtree& tr,vtree::iterator it) const {
    typedef std::multiset<pre_it, util::lexicographic_subtree_order<vertex>
      > subtree_quotient;
    typedef subtree_quotient::iterator subtree_quotient_it;
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      if(is_builtin(*b1)) {
	builtin b1_builtin = get_builtin(*b1);
	if(b1_builtin==*b2) {
	  if(b1_builtin==id::div) {
	    assert(tr.number_of_children(b1)==2);
	    assert(tr.number_of_children(b2)==2);
	    pre_it b1_num = tr.child(b1, 0);
	    pre_it b1_denum = tr.child(b1, 1);
	    pre_it b2_num = tr.child(b2, 0);
	    pre_it b2_denum = tr.child(b2, 1);
	    if(tr.equal_subtree(b1_num, b2_num)) {
	      *it = id::div;
	      pre_it new_if = tr.prepend_child(it, id::contin_if);
	      tr.reparent(new_if, ++(it.begin()), it.end());
	      tr.move_before(new_if, b1_num);
	      tr.erase(b2_num);
	      tr.erase(tr.flatten(b1));
	      tr.erase(tr.flatten(b2));
	    }
	    else if(tr.equal_subtree(b1_denum, b2_denum)) {
	      *it = id::div;
	      pre_it new_if = tr.prepend_child(it, id::contin_if);
	      tr.reparent(new_if, ++(it.begin()), it.end());
	      tr.move_after(new_if, b1_denum);
	      tr.erase(b2_denum);
	      tr.erase(tr.flatten(b1));
	      tr.erase(tr.flatten(b2));
	    }
	  }
	  if(b1_builtin==id::plus || b1_builtin==id::times) {
	    subtree_quotient sq;
	    pre_it new_if;
	    //fill sq
	    for(sib_it sib = b1.begin(); sib != b1.end(); ++sib)
	      sq.insert(pre_it(sib));
	    //move redundant arguments of op out of contin_if
	    subtree_quotient_it sq_it;
	    for(sib_it sib = b2.begin(); sib != b2.end();) {
	      sq_it = sq.find(pre_it(sib));
	      if(sq_it!=sq.end()) {
		if(!tr.is_valid(new_if)) {
		  *it = b1_builtin;
		  new_if = tr.prepend_child(it, id::contin_if);
		  tr.reparent(new_if, ++(it.begin()), it.end());
		}
		tr.move_after(new_if, *sq_it);
		sib = tr.erase(sib);
		sq.erase(sq_it);
	      }
	      else ++sib;
	    }
	    //check if old op is empty or contin_if is empty
	    if(b1.has_one_child())
	      tr.erase(tr.flatten(b1));
	    else if(b1.is_childless())
	      *b1 = (b1_builtin==id::plus?0.0:1.0);
	    if(b2.has_one_child())
	      tr.erase(tr.flatten(b2));
	    else if(b2.is_childless())
	      *b2 = (b1_builtin==id::plus?0.0:1.0);
	  }
	}
      }
    }
  }

  //contin_if(x exp1[x] exp2[x]) -> if(x exp1[true] exp2[false])
  void reduce_contin_if_substitute_cond::operator()(vtree& tr,
						    vtree::iterator it) const {
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it cond = tr.child(it, 0);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      if(!may_have_side_effects(cond)) {
	if(*cond!=id::logical_true && *cond!=id::logical_false) {
	  pre_it cond_in_b1 =
	    tr.find_subtree(cond, pre_it(b1.begin()), b2);
	  while(tr.is_valid(cond_in_b1)) {
	    *cond_in_b1 = id::logical_true;
	    tr.erase_children(cond_in_b1);
	    cond_in_b1 =
	      tr.find_subtree(cond, ++cond_in_b1, b2);
	  }
	  pre_it cond_in_b2 =
	    tr.find_subtree(cond, pre_it(b2.begin()), pre_it(b2.end()));
	  while(tr.is_valid(cond_in_b2)) {
	    *cond_in_b2 = id::logical_false;
	    tr.erase_children(cond_in_b2);
	    cond_in_b2 =
	      tr.find_subtree(cond, ++cond_in_b2, pre_it(b2.end()));
	  }
	}
      }
    }
  }

  //0<x+c1 and 0<x+c2 -> 0<x+min(c1, c2)
  //or more generally
  //0<sum x_i +c1 and 0<sum x_i +c2 and 0<sum x_i +c3...
  //->sum x_i + min(c1, c2, c3, ...)
  //
  //and
  //
  //0<x+c1 or 0<x+c2 -> 0<x+max(c1, c2)
  //or more generally
  //0<sum x_i +c1 and 0<sum x_i +c2 and 0<sum x_i +c3...
  //->sum x_i + max(c1, c2, c3, ...)
  void reduce_junction_gt_zero_sum_constant::operator()(vtree& tr,
							vtree::iterator it) const {
    //the first pre_it is the root of the subtree s obtained from
    //0<s+c
    //such that s is a sum of non-constant term and c is a constant
    //the second argument, contin_t, is c
    typedef std::map<pre_it, contin_t,
      util::lexicographic_subtree_order<vertex>
      > subtree_quotient;
    typedef subtree_quotient::iterator subtree_quotient_it;
    if(*it==id::logical_and || *it==id::logical_or) {
      subtree_quotient sq;
      for(sib_it sib = it.begin(); sib != it.end();) {
	if(*sib==id::greater_than_zero) {
	  assert(sib.has_one_child());
	  pre_it sib_child = sib.begin();
	  if(*sib_child==id::plus) {
	    contin_t c = 0.0;
	    for(sib_it plus_sib = sib_child.begin(); plus_sib != sib_child.end();)
	      if(is_contin(*plus_sib)) {
		c = c + get_contin(*plus_sib);
		plus_sib = tr.erase(plus_sib);
	      }
	      else ++plus_sib;
	    if(sib_child.is_childless()) {
	      *sib_child = c;
	      ++sib;
	    }
	    else {
	      subtree_quotient_it sqi = sq.find(pre_it(sib_child));
	      if(sqi==sq.end()) {
		std::pair<pre_it, contin_t> p(pre_it(sib_child), c);
		sq.insert(p);
		++sib;
	      }
	      else {
		contin_t& c_ref = sqi->second;
		c_ref = (*it==id::logical_and?
			 std::min(c, c_ref) : std::max(c, c_ref));
		sib = tr.erase(sib);
	      }
	    }
	  }
	  else ++sib;
	}
	else ++sib;
      }
      //add the right constant [that is min/max(c1, c2, ....)] to each
      //reduced inequality
      for(subtree_quotient_it i = sq.begin(); i != sq.end(); ++i) {
	contin_t val = i->second;
	if(val != 0.0)
	  tr.append_child(i->first, i->second);
      }
    }
  }

  //implies
  //takes a tree and 2 iterators and tells if it1 implies it2
  //for instance x implies x, because x==x
  //if x!=y for instance x==2<z implies y==3<z
  bool reduce_from_assumptions::implies(const vtree& tr,
					vtree::iterator it1,
					vtree::iterator it2) const {
    //if there are syntactly equal clearly it1 implies it2
    if(tr.equal_subtree(it1, it2))
      return true;
    //otherwise deal with inequalities
    else {
      pre_it ine1_child; //child of inequality ot it1
      bool is_strict1;
      pre_it ine2_child; //child of inequality of it2
      bool is_strict2;
      //determine if ine1_child is a strict inequality
      if(*it1==id::greater_than_zero) {
	assert(it1.has_one_child());
	ine1_child = it1.begin();
	is_strict1 = true;
      }
      //determine if ine1_child is a non-strict inequality
      else if(*it1==id::logical_not) {
	assert(it1.has_one_child());
	pre_it it1_child = it1.begin();
	if(*it1_child==id::greater_than_zero) {
	  assert(it1_child.has_one_child());
	  ine1_child = it1_child.begin();
	  is_strict1 = false;
	}
      }
      //determine if ine2_child is a strict inequality
      if(*it2==id::greater_than_zero) {
	assert(it2.has_one_child());
	ine2_child = it2.begin();
	is_strict2 = true;
      }
      //determine if ine2_child is a non-strict inequality
      else if(*it2==id::logical_not) {
	assert(it2.has_one_child());
	pre_it it2_child = it2.begin();
	if(*it2_child==id::greater_than_zero) {
	  assert(it2_child.has_one_child());
	  ine2_child = it2_child.begin();
	  is_strict2 = false;
	}
      }
      //check is it1 and it2 are both inequalities
      if(tr.is_valid(ine1_child) && tr.is_valid(ine2_child)) {
	contin_t c1 = 0.0;
	contin_t c2 = 0.0;
	//tr1 will contain the non constant terms of inequality 1
	vtree tr1 = tr.subtree(ine1_child, tr.next_sibling(ine1_child));
	pre_it tr1_it = tr1.begin();
	//tr2 will contain the non constant terms of inequality 2
	vtree tr2 = tr.subtree(ine2_child, tr.next_sibling(ine2_child));
	pre_it tr2_it = tr2.begin();

	//the idea is that all constant under + will be take off from and 
	//tr1 and tr2 then tr1 and tr2 will be compared if there are equal
	//(but previously non strict equality will be multiplied
	//by -1 because non strict inequality is represented using not(0<...) )
	//because of this possible multiplication tr1 and tr2 will be first
	//reduced.

	//take off constant terms from tr1 and put in c1
	//if non strict c1 is inverted because non strict are represented
	//with not(0<...)
	if(*tr1_it==id::plus) {
	  for(sib_it sib1 = tr1_it.begin(); sib1 != tr1_it.end();) {
	    if(is_contin(*sib1)) {
	      c1 = c1 + (is_strict1? 1.0 : -1.0) * get_contin(*sib1);
	      sib1 = tr1.erase(sib1);
	    }
	    else ++sib1;
	  }
	}
	//check if tr1 under the form -1*(x1+..+xn+c) and take off c
	else if(*tr1_it==id::times && tr1_it.number_of_children()==2) {
	  pre_it child0 = tr1.child(pre_it(tr1_it), 0);
	  pre_it child1 = tr1.child(pre_it(tr1_it), 1);
	  //check if -1 in first child
	  if(is_contin(*child0) && get_contin(*child0)==-1.0) {
	    //take off c
	    if(*child1==id::plus) {
	      for(sib_it s1p = child1.begin(); s1p != child1.end();)
		if(is_contin(*s1p)) {
		  //this time c is not inverted if non strict because already
		  //inverted with -1*...
		  c1 = c1 + (is_strict1? -1.0 : 1.0) * get_contin(*s1p);
		  s1p = tr1.erase(s1p);
		}
	        else ++s1p;
	    }
	  }
	  //check if -1 in the second child
	  else if(is_contin(*child1) && get_contin(*child1)==-1.0) {
	    //take off c
	    if(*child0==id::plus) {
	      for(sib_it s1p = child0.begin(); s1p != child0.end();)
		if(is_contin(*s1p)) {
		  //this time c is not inverted if non strict because already
		  //inverted with -1*...
		  c1 = c1 + (is_strict1? -1.0 : 1.0) * get_contin(*s1p);
		  s1p = tr1.erase(s1p);
		}
	        else ++s1p;
	    }
	  }
	}
	if(!is_strict1) //multiply tr1 by -1
	  tr1.append_child(tr1.wrap(tr1_it, id::times), -1.0);
	
	//take off constant terms from tr2 and put in c2
	//if non strict c2 is inverted because non strict are represented
	//with not(0<...)
	if(*tr2_it==id::plus) {
	  for(sib_it sib2 = tr2_it.begin(); sib2 != tr2_it.end();) {
	    if(is_contin(*sib2)) {
	      c2 = c2 + (is_strict2? 1.0 : -1.0) * get_contin(*sib2);
	      sib2 = tr2.erase(sib2);
	    }
	    else ++sib2;
	  }
	}
	//check if tr2 under the form -1*(x1+..+xn+c) and take off c
	else if(*tr2_it==id::times && tr2_it.number_of_children()==2) {
	  pre_it child0 = tr2.child(pre_it(tr2_it), 0);
	  pre_it child1 = tr2.child(pre_it(tr2_it), 1);
	  //check if -1 in first child
	  if(is_contin(*child0) && get_contin(*child0)==-1.0) {
	    //take off c
	    if(*child1==id::plus) {
	      for(sib_it s2p = child1.begin(); s2p != child1.end();)
		if(is_contin(*s2p)) {
		  //this time c is not inverted if non strict because already
		  //inverted with -1*...
		  c2 = c2 + (is_strict2? -1.0 : 1.0) * get_contin(*s2p);
		  s2p = tr2.erase(s2p);
		}
	        else ++s2p;
	    }
	  }
	  //check if -1 in the second child
	  else if(is_contin(*child1) && get_contin(*child1)==-1.0) {
	    //take off c
	    if(*child0==id::plus) {
	      for(sib_it s2p = child0.begin(); s2p != child0.end();)
		if(is_contin(*s2p)) {
		  //this time c is not inverted if non strict because already
		  //inverted with -1*...
		  c2 = c2 + (is_strict2? -1.0 : 1.0) * get_contin(*s2p);
		  s2p = tr2.erase(s2p);
		}
	        else ++s2p;
	    }
	  }
	}
	if(!is_strict2) //multiply tr1 by -1
	  tr1.append_child(tr2.wrap(tr2_it, id::times), -1.0);

	//reduce tr1 and tr2 and compare is equal
	(*_reduction)(tr1);
	(*_reduction)(tr2);

	if(tr1==tr2) {
	  if(!is_strict1 && is_strict2)
	    return c1 < c2;
	  else return c1 <= c2;
	}
	else return false;
      }
      //no equality, no inequalities
      else return false;
    }
  }
  
  //look up the assumptions and replace by true is present or false
  //if not(assum) present
  void reduce_from_assumptions::operator()(vtree& tr,vtree::iterator it) const {
    if((is_boolean_output(*it)
	&& *it!=id::logical_true && *it!=id::logical_false)
       || is_argument(*it)) {
      sib_it assum = tr.next_sibling(tr.begin());
      if(tr.is_valid(assum)) {
	for(; assum != tr.end(); ++assum) {
	  if(implies(tr, pre_it(assum), it)) {
	    *it = id::logical_true;
	    tr.erase_children(it);
	    return;
	  }
	  else if(*assum==id::logical_not) {
	    assert(assum.has_one_child());
	    if(implies(tr, it, pre_it(assum.begin()))) {
	      *it = id::logical_false;
	      tr.erase_children(it);
	      return;
	    }
	  }
	  else if(is_argument(*assum) && is_argument(*it)) {
	    argument assum_arg = get_argument(*assum);
	    assum_arg.negate();
	    //don't need to treat the case when same idx because
	    //already treated in the beginning of the function
	    if(assum_arg.idx == get_argument(*it).idx) {
	      *it = id::logical_false;
	      tr.erase_children(it);
	      return;
	    }
	  }
	}
      }
    }
  }

  //if(x y z) -> if(NOT(x) z y)  if |NOT(x)|<|x|
  void reduce_contin_if_not::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::contin_if) {
      assert(tr.number_of_children(it)==3);
      pre_it cond = tr.child(it, 0);
      pre_it b1 = tr.child(it, 1);
      pre_it b2 = tr.child(it, 2);
      vtree cond_tr = tr.subtree(sib_it(cond), sib_it(b1));
      //copy old assumptions, begin
      sib_it bna = cond_tr.begin(); //before new assumption
      for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	bna = cond_tr.insert_subtree_after(bna, a);
      //copy old assumptions, end
      pre_it cond_tr_it = cond_tr.begin();
      pre_it new_cond_it = cond_tr.prepend_child(cond_tr_it, *cond_tr_it);
      sib_it first_cond_child = cond_tr.next_sibling(sib_it(new_cond_it));
      *cond_tr_it = id::logical_not;
      cond_tr.reparent(new_cond_it, first_cond_child, cond_tr_it.end());
      (*_reduction)(cond_tr);
      (*_reduction)(tr, cond);
      if(cond_tr.size() < tr.subtree_size(cond)) {
	tr.replace(sib_it(cond), sib_it(b1),
		   sib_it(cond_tr.begin()), sib_it(cond_tr.end()));
	tr.swap(b1, b2);
      }
    }
  }

  //0<sum x_i -> true    if 0<x_i -> true forall i
  //0<sum x_i -> false   if 0<x_i -> false forall i
  void reduce_gt_zero_sum::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::plus) {
	bool res;
	for(sib_it sib = it_child.begin(); sib != it_child.end(); ++sib) {
	  vtree sib_tr = tr.subtree(sib, tr.next_sibling(sib));
	  //copy old assumptions, begin
	  sib_it bna = sib_tr.begin(); //before new assumption
	  for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	    bna = sib_tr.insert_subtree_after(bna, a);
	  //copy old assumptions, end
	  sib_tr.wrap(sib_tr.begin(), id::greater_than_zero);
	  (*_reduction)(sib_tr);
	  pre_it sib_tr_root = sib_tr.begin();
	  if(is_boolean(*sib_tr_root)) {
	    if(tr.sibling_index(sib)==0)
	      res = vertex_to_bool(*sib_tr_root);
	    else if(res != vertex_to_bool(*sib_tr_root))
	      return;
	  }
	  else return;
	}
	*it = bool_to_vertex(res);
	tr.erase_children(it);
      }
    }
  }

  //0<prod x_i -> 0<prod x_i    with x_i=1 if 0<x_i -> true
  //                            with x_i=-1 if 0<-x_i -> true
  //0<prod x_i -> false     if exist x_i==0, 
  //                           that is 0<x_i -> false and 0<-1*x_i -> false
  void reduce_gt_zero_prod::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::times) {
	for(sib_it sib = it_child.begin(); sib != it_child.end();) {
	  vtree sib_tr = tr.subtree(sib, tr.next_sibling(sib));
	  //copy old assumptions, begin
	  sib_it bna = sib_tr.begin(); //before new assumption
	  for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	    bna = sib_tr.insert_subtree_after(bna, a);
	  //copy old assumptions, end
	  sib_tr.wrap(sib_tr.begin(), id::greater_than_zero);
	  //reduce 0<x_i
	  vtree tr1 = sib_tr;
	  (*_complete_reduction)(tr1);
	  pre_it tr1_root = tr1.begin();
	  if(*tr1_root==id::logical_true)
	    sib = tr.erase(sib);
	  else {
	    //reduce 0<-1*x_i
	    vtree tr2 = sib_tr;
	    pre_it x_i = tr2.begin().begin();
	    tr2.append_child(tr2.wrap(x_i, id::times), -1.0);
	    //reduction without reduce_gt_zero_prod is used to avoid infinite
	    //recursion
	    (*_reduction_without_itself)(tr2);
	    pre_it tr2_root = tr2.begin();
	    if(*tr2_root==id::logical_true) {
	      *sib = -1.0;
	      tr.erase_children(sib);
	      ++sib;
	    }
	    //test if x_i==0.0
	    else if(*tr1_root==id::logical_false
		    && *tr2_root==id::logical_false) {
	      *it = id::logical_false;
	      tr.erase_children(it);
	      return;
	    }
	    else ++sib;
	  }
	}
	if(it_child.is_childless()) {
	  tr.erase_children(it);
	  *it = id::logical_true;
	}
      }
    }
  }

  //0<x/y -> 0<x    if 0<y -> true
  //0<x/y -> 0<y    if 0<x -> true
  //0<x/y -> 0<-1*y if 0<-1*x -> true
  //0<x/y -> 0<-1*x if 0<-1*y -> true
  //0<x/y -> false  if x==0, that is not(0<x) -> true and not(0<-x) -> true
  void reduce_gt_zero_div::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::div) {
	assert(it_child.number_of_children()==2);
	pre_it num = tr.child(it_child, 0);
	pre_it denom = tr.child(it_child, 1);
	//case of the numerator
	vtree num_tr = tr.subtree(sib_it(num), sib_it(denom));
	//copy old assumptions, begin
	sib_it bna = num_tr.begin(); //before new assumption
	for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	  bna = num_tr.insert_subtree_after(bna, a);
	//copy old assumptions, end
	num_tr.wrap(num_tr.begin(), id::greater_than_zero);
	//reduce 0<x
	vtree num_tr1 = num_tr;
	(*_reduction)(num_tr1);
	pre_it num_tr1_root = num_tr1.begin();
	if(*num_tr1_root==id::logical_true) {
	  tr.erase(num);
	  tr.erase(tr.flatten(it_child));
	}
	else {
	  //reduce 0<-1*x
	  vtree num_tr2 = num_tr;
	  num_tr2.append_child(tr.wrap(num_tr2.begin().begin(), id::times),
			       -1.0);
	  (*_reduction)(num_tr2);
	  pre_it num_tr2_root = num_tr2.begin();
	  if(*num_tr2_root==id::logical_true) {
	    tr.erase_children(num);
	    *it_child = id::times;
	    *num = -1.0;
	  }
	  //test if x==0
	  else if(*num_tr1_root==id::logical_false
		  && *num_tr2_root==id::logical_false) {
	    *it = id::logical_false;
	    tr.erase_children(it);
	  }
	  //case of denominateur
	  else {
	    vtree denom_tr = tr.subtree(sib_it(denom), tr.next_sibling(denom));
	    //copy old assumptions, begin
	    pre_it bna = denom_tr.begin(); //before new assumption
	    for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	      bna = denom_tr.insert_subtree_after(bna, a);
	    //copy old assumptions, end
	    denom_tr.wrap(denom_tr.begin(), id::greater_than_zero);
	    //reduce 0<y
	    vtree denom_tr1 = denom_tr;
	    (*_reduction)(denom_tr1);
	    pre_it denom_tr1_root = denom_tr1.begin();
	    if(*denom_tr1_root==id::logical_true) {
	      tr.erase(denom);
	      tr.erase(tr.flatten(it_child));
	    }
	    else {
	      //reduce 0<-1*y
	      vtree denom_tr2 = denom_tr;
	      denom_tr2.append_child(tr.wrap(denom_tr2.begin().begin(),
					     id::times), -1.0);
	      (*_reduction)(denom_tr2);
	      pre_it denom_tr2_root = denom_tr2.begin();
	      if(*denom_tr2_root==id::logical_true) {
		tr.erase_children(denom);
		*it_child = id::times;
		*denom = -1.0;
	      }
	    }
	  }
	}
      }
    }
  }

  //0<x+sin(y) -> true  if 0<x-1 -> true
  //0<x+sin(y) -> false if 0<x+1 -> false
  void reduce_gt_zero_sum_sin::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::plus && tr.is_valid(it_child.find_child(id::sin))) {
	vtree copy_tr = tr.subtree(sib_it(it), tr.next_sibling(sib_it(it)));
	//copy old assumptions, begin
	sib_it bna = copy_tr.begin(); //before new assumption
	for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	  bna = copy_tr.insert_subtree_after(bna, a);
	//copy old assumptions, end
	pre_it copy_plus = copy_tr.begin().begin();
	int number_of_sin = 0;
	for(sib_it sib = copy_plus.begin(); sib != copy_plus.end();) {
	  if(*sib==id::sin) {
	    number_of_sin++;
	    sib = copy_tr.erase(sib);
	  }
	  else ++sib;
	}
	if(copy_plus.is_childless())
	  return;
	//check if 0<x+y..-number_of_sin -> true
	vtree copy1_tr = copy_tr;
	pre_it copy1_plus = copy1_tr.begin().begin();
	copy1_tr.append_child(copy1_plus, (contin_t)-number_of_sin);
	(*_reduction)(copy1_tr);
	if(*copy1_tr.begin()==id::logical_true) {
	  *it = id::logical_true;
	  tr.erase_children(it);
	}
	else {
	  //check if 0<x+y..+number_of_sin -> false
	  vtree copy2_tr = copy_tr;
	  pre_it copy2_plus = copy2_tr.begin().begin();
	  copy2_tr.append_child(copy2_plus, (contin_t)number_of_sin);
	  (*_reduction)(copy2_tr);
	  if(*copy2_tr.begin()==id::logical_false) {
	    *it = id::logical_false;
	    tr.erase_children(it);
	  }
	}
      }
    }
  }

  //0<sin(y) -> true  if 0<y -> true and 0<pi-y -> true
  //0<sin(y) -> false if 0<y -> false and 0<-(y+pi) -> false
  void reduce_gt_zero_sin::operator()(vtree& tr,vtree::iterator it) const {
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::sin) {
	assert(it_child.has_one_child());
	pre_it y = it_child.begin();
	vtree copy_tr = tr.subtree(sib_it(it), tr.next_sibling(sib_it(it)));
	//copy old assumptions, begin
	sib_it bna = copy_tr.begin(); //before new assumption
	for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	  bna = copy_tr.insert_subtree_after(bna, a);
	//copy old assumptions, end
	//turn copy_tr into 0<y and check if 0<y
	copy_tr.erase(copy_tr.flatten(copy_tr.begin().begin()));
	(*_reduction)(copy_tr);
	if(*copy_tr.begin()==id::logical_true) {
	  //check if 0<pi-y -> true
	  vtree copy1_tr = tr.subtree(sib_it(it), tr.next_sibling(sib_it(it)));
	  //copy old assumptions, begin
	  sib_it bna = copy1_tr.begin(); //before new assumption
	  for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	    bna = copy1_tr.insert_subtree_after(bna, a);
	  //copy old assumptions, end
	  pre_it sin1_it = copy1_tr.begin().begin();
	  *sin1_it = id::times;
	  copy1_tr.append_child(sin1_it, -1.0);
	  pre_it new_plus = copy1_tr.wrap(sin1_it, id::plus);
	  copy1_tr.append_child(new_plus, (contin_t)PI);
	  (*_reduction)(copy1_tr);
	  if(*copy1_tr.begin()==id::logical_true) {
	    *it = id::logical_true;
	    tr.erase_children(it);
	  }
	}
	else if(*copy_tr.begin()==id::logical_false) {
	  //check if 0<-(y+pi) -> false //TODO
	  vtree copy2_tr = tr.subtree(sib_it(it), tr.next_sibling(sib_it(it)));
	  //copy old assumptions, begin
	  sib_it bna = copy2_tr.begin(); //before new assumption
	  for(sib_it a = tr.next_sibling(tr.begin()); a != tr.end(); ++a)
	    bna = copy2_tr.insert_subtree_after(bna, a);
	  //copy old assumptions, end
	  pre_it sin2_it = copy2_tr.begin().begin();
	  *sin2_it = id::plus;
	  copy2_tr.append_child(sin2_it, (contin_t)PI);
	  copy2_tr.append_child(copy2_tr.wrap(sin2_it, id::times), -1.0);
	  (*_reduction)(copy2_tr);
	  if(*copy2_tr.begin()==id::logical_false) {
	    *it = id::logical_false;
	    tr.erase_children(it);
	  }
	}
      }
    }
  }

  
  contin_t reduce_inequality_from_assumptions::splitCoefExpression(vtree& tr, vtree::iterator& it) const {
    contin_t coef = 1.0;
    if(*it==id::times) {
      for(sib_it sib = it.begin(); sib != it.end();) {
	if(is_contin(*sib)) {
	  coef = coef * get_contin(*sib);
	  sib = tr.erase(sib);
	}
	else ++sib;
      }
      if(it.has_one_child())
	it = tr.erase(tr.flatten(it));
      else if(it.is_childless()) {
	tr.erase(it);
	it = tr.end(); //to be recognized as not valid
      }
    }
    else if(is_contin(*it)) {
      coef = get_contin(*it);
      tr.erase(it);
      it = tr.end(); //to be recognized as not valid
    }
    return coef;
  }

  bool reduce_inequality_from_assumptions::gaussJordanElimination(double_matrix& dm, double eps) const {
    unsigned int h = dm.size(); //number of rows
    assert(h>0);
    unsigned int w = dm[h-1].size(); //number of columns
    for(unsigned int y = 0; y < h; y++) {
      unsigned int maxrow = y;
      //find max pivot
      for(unsigned int y2 = y+1; y2 < h; y2++)
	if(std::fabs(dm[y2][y]) > std::fabs(dm[maxrow][y]))
	  maxrow = y2;
      std::vector<double> tmp_row = dm[y];
      dm[y] = dm[maxrow];
      dm[maxrow] = tmp_row;
      //check if singular
      if(std::fabs(dm[y][y]) <= eps)
	return false;
      //eliminate column y
      for(unsigned int y2 = y+1; y2 < h; y2++) {
	double c = dm[y2][y]/dm[y][y];
	for(unsigned int x = y; x < w; x++)
	  dm[y2][x] -= c*dm[y][x];
      }
    }
    //backsubstitute
    for(int y = h-1; y >= 0; y--) {
      double c = dm[y][y];
      for(int y2 = 0; y2 < y; y2++)
	for(int x = w-1; x > y-1; x--)
	  dm[y2][x] -= dm[y][x]*dm[y2][y]/c;
      dm[y][y] /= c;
      for(unsigned x = (h>w?w:h); x < (h>w?h:w); x++)
	dm[y][x] /= c;
    }
    return true;
  }

  void reduce_inequality_from_assumptions::gaussianElimination(double_matrix& dm) const {
    unsigned int m = dm.size(); //number of rows
    assert(m>0);
    unsigned int n = dm[m-1].size(); //number of columns
    for(unsigned int i = 0; i < m; i++) {
      for(unsigned int j = 0; j < n; j++) {
	unsigned int maxi = i;
	for(unsigned int k = i+1; k < m; k++) {
	  if(std::fabs(dm[k][j])>fabs(dm[maxi][j]))
	    maxi = k;
	}
	if(dm[maxi][j] != 0.0) {
	  std::vector<double> tmp_row = dm[i];
	  dm[i] = dm[maxi];
	  dm[maxi] = tmp_row;
	  double c = dm[i][j];
	  assert(c != 0.0);
	  for(unsigned int k = 0; k < n; k++)
	    dm[i][k] = dm[i][k] / c;
	  for(unsigned int u = i+1; u < m; u++) {
	    double d = dm[u][j];
	    for(unsigned int k = 0; k < n; k++)
	      dm[u][k] = dm[u][k] - d*dm[i][k];
	  }
	}
      }
    }
  }

  bool reduce_inequality_from_assumptions::gaussianElimination2(double_matrix& dm, double eps) const {
    int m = dm.size(); //number of rows
    assert(m>0);
    int n = dm[m-1].size(); //number of columns
    assert(n>0);
    for(int i = 0; i < std::min(m, n-1); i++) {
      //swap below until find non-zero
      int ipivot = -1;
      for(int i2  = i; i2 < m; i2++)
	if(std::fabs(dm[i2][i]) > eps)
	  ipivot = i2;
      //swap row i and row
      if(ipivot==-1)
	return false;
      std::vector<contin_t> tmp = dm[i];
      dm[i] = dm[ipivot];
      dm[ipivot] = tmp;
      //divide pivot by itself to be 1
      contin_t c = dm[i][i];
      for(int j = 0; j < n; j++)
	dm[i][j] /= c;
      //reduce to 0 all column below
      for(int i2 = i+1; i2 < m; i2++) {
	contin_t c = dm[i2][i];
	for(int j = 0; j < n; j++)
	  dm[i2][j] -= c*dm[i][j];
      }
    }
    return true;
  }

  bool reduce_inequality_from_assumptions::backSubstitution(double_matrix& dm, std::vector<double>& solution, double eps) const {
    assert(solution.empty());
    int m = dm.size(); //number of rows
    assert(m>0);
    int n = dm[m-1].size(); //number of columns
    assert(n>0);
    solution.resize(n-1, 0.0);
    bool allNull = true;
    for(int i = m-1; i >= 0; i--) {
      //check if there is contradiction in the equationnal system
      for(int j = 0; j < n-1; j++)
	if(std::fabs(dm[i][j]) > eps) 
	  allNull = false;
      if(allNull && std::fabs(dm[i][n-1]) > eps) //something in contradiction
	return false;
      for(int j = 0; j < n-1; j++) {	
	if(dm[i][j]==1.0) {
	  solution[j] = dm[i][n-1];
	  for(int i2 = i-1; i2 >= 0; i2--) {
	    dm[i2][n-1] -= solution[j]*dm[i2][j];
	    dm[i2][j] = 0;
	  }
	}
      }
    }
    return true;
  }
							    

  void reduce_inequality_from_assumptions::operator()
    (vtree& tr,vtree::iterator it) const {
    //std::cout << "LINEAR BEFORE : " << tr << " IT : " << *it << std::endl;
    //associate a subtree and row in the matrix to solve
    //the size of a row is |assumptions|+1
    //the last column of the matrix contains the coef associated with its
    //expression
    typedef std::map<pre_it, std::vector<contin_t>,
      util::lexicographic_subtree_order<vertex> > subtree_row;
    typedef subtree_row::iterator subtree_row_it;
    //true if the main term is strictly positive false otherwise
    bool termStrictPositive;
    //we make a copy of tr because the reductor is going
    //manipulate/change it to find expression/variables.
    vtree ctr = tr;
    //main term to find the linear decomposition
    pre_it term_it;
    //check if strict positive inequality or negative or null
    if(*it==id::greater_than_zero) {
      assert(it.has_one_child());
      term_it = ctr.begin().begin();
      termStrictPositive = true;
    }
    else if(*it==id::logical_not) {
      assert(it.has_one_child());
      pre_it it_child = it.begin();
      if(*it_child==id::greater_than_zero) {
	assert(it_child.has_one_child());
	term_it = ctr.begin().begin().begin();
	termStrictPositive = false;
      }
    }
    //check if the main term is an inequality (positive or negative)
    if(ctr.is_valid(term_it)) {
      sib_it assum = ctr.next_sibling(ctr.begin());
      if(ctr.is_valid(assum)) {
	//select the set of inequality assumptions
	//if this set is empty then stop the reduction
	std::vector<pre_it> assumptions;
	//define a boolean vector such that
	//has ther size of number of assumptions + 1 (for the main term)
	//(the main term is at the last position
	//if an entry if true it means that the assumption (or the last term)
	//is strictly positive
	//otherwise it means that it is negative or null
	std::vector<bool> isStrictPositives;
	for(; assum != ctr.end(); ++assum) {
	  if(*assum==id::greater_than_zero) {
	    assumptions.push_back(pre_it(assum));
	    isStrictPositives.push_back(true);
	  }
	  else if(*assum==id::logical_not) {
	    assert(assum.has_one_child());
	    pre_it assum_child = assum.begin();
	    if(*assum_child==id::greater_than_zero) {
	      assumptions.push_back(assum_child);
	      isStrictPositives.push_back(false);
	    }
	  }
	}
	//the last entry of true because a constant is consider
	//strictly positive of a term 1
	isStrictPositives.push_back(true);
	if(!assumptions.empty()) {
	  subtree_row sr;
	  //coefs alone without term
	  std::vector<contin_t> free_coefs(assumptions.size()+2, 0.0);
	  //for the constant term in the linear combination
	  free_coefs[assumptions.size()] = 1.0;
	  //this reduction only handle linear combination
	  //that's why the case when term_it is * is not treated
	  //if the main term is c*+(...) with c positive
	  //then c would be eliminated by previous reduction rules
	  //however if c is negative it is then replaced by -1
	  //so the only case that is treated when term_it is * is
	  //-1*x
	  contin_t main_coef = 1.0; //to keep in mind that they were -1*term_it
	  if(*term_it==id::times && term_it.number_of_children()==2) {
	    pre_it child0 = ctr.child(pre_it(term_it), 0);
	    pre_it child1 = ctr.child(pre_it(term_it), 1);
	    if(is_contin(*child0) && get_contin(*child0)==-1.0) {
	      term_it = child1;
	      main_coef = -1.0;
	    }
	    else if(is_contin(*child1) && get_contin(*child1)==-1.0) {
	      term_it = child0;
	      main_coef = -1.0;
	    }
	  }
	  //std::cout << "MAIN COEF : " << main_coef << std::endl;
	  if(*term_it==id::plus) {
	    for(sib_it sib = term_it.begin(); sib != term_it.end();) {
	      //identify paires coef*expression
	      pre_it expression = sib;
	      sib_it next_sib = ctr.next_sibling(sib);
	      contin_t coef = splitCoefExpression(ctr, expression) * main_coef;
	      sib = next_sib;
	      //add or uptdate (expression, coef) in the map
	      if(ctr.is_valid(expression)) {
		subtree_row_it sr_it = sr.find(expression);
		if(sr_it==sr.end()) {
		  std::vector<contin_t> v(assumptions.size()+1, 0.0);
		  v.push_back(coef);
		  std::pair<pre_it, std::vector<contin_t> > p(expression, v);
		  sr.insert(p);
		}
		else { //update coef
		  std::vector<contin_t>& v_sr = sr_it->second;
		  v_sr[assumptions.size()+1] = v_sr[assumptions.size()+1]+coef;
		}
	      }
	      //add coef to the last entry of free_coefs
	      else free_coefs[assumptions.size()+1] 
		     = free_coefs[assumptions.size()+1] + coef;
	      //print map
	      /*for(subtree_row_it sr_it = sr.begin(); sr_it != sr.end(); ++sr_it) {
		std::cout << "MAIN VAR : " << *sr_it->first << " ";
		for(std::vector<double>::iterator i = (sr_it->second).begin();
		    i != (sr_it->second).end(); ++i) {
		  std::cout << *i << " ";
		}
		std::cout << std::endl;
		}*/
	    }
	  }
	  //main term is not a sum
	  else {
	    //identify paires coef*expression
	    pre_it expression = term_it;
	    contin_t coef = splitCoefExpression(ctr, expression) * main_coef;
	    if(ctr.is_valid(expression)) {
	      std::vector<contin_t> v(assumptions.size()+1, 0.0);
	      v.push_back(coef);
	      std::pair<pre_it, std::vector<contin_t> > p(expression, v);
	      sr.insert(p);
	    }
	    //add coef to the last entry of free_coefs
	    else free_coefs[assumptions.size()+1] 
		   = free_coefs[assumptions.size()+1] + coef;
	  }
	  //complete sr and free_coefs with the assumptions
	  //either by adding new expressions not found in ctr and coef 0.0
	  //or by filling the matrix
	  for(unsigned int i = 0; i < assumptions.size(); i++) {
	    pre_it assum = assumptions[i];
	    assert(*assum==id::greater_than_zero && assum.has_one_child());
	    pre_it assum_child = assum.begin();
	    //same remark as previously about -1*x
	    contin_t assum_coef = 1.0;
	    if(*assum_child==id::times && assum_child.number_of_children()==2){
	      pre_it child0 = ctr.child(pre_it(assum_child), 0);
	      pre_it child1 = ctr.child(pre_it(assum_child), 1);
	      if(is_contin(*child0) && get_contin(*child0)==-1.0) {
		assum_child = child1;
		assum_coef = -1.0;
	      }
	      else if(is_contin(*child1) && get_contin(*child1)==-1.0) {
		assum_child = child0;
		assum_coef = -1.0;
	      }
	    }
	    if(*assum_child==id::plus) {
	      for(sib_it pl_sib = assum_child.begin();
		  pl_sib != assum_child.end();) {
		//identify paires coef*expression
		pre_it expression = pl_sib;
		sib_it next_pl_sib = ctr.next_sibling(pl_sib);
		contin_t coef = 
		  splitCoefExpression(ctr, expression) * assum_coef;
		pl_sib = next_pl_sib;
		//add or uptdate (expression, rows) in the map
		if(ctr.is_valid(expression)) {
		  subtree_row_it sr_it = sr.find(expression);
		  if(sr_it==sr.end()) {
		    std::vector<contin_t> v(assumptions.size()+2, 0.0);
		    v[i] = coef;
		    std::pair<pre_it, std::vector<contin_t> >
		      p(expression, v);
		    sr.insert(p);
		  }
		  else { //update coef
		    std::vector<contin_t>& v_sr = sr_it->second;
		    v_sr[i] = v_sr[i] + coef;
		  }
		}
		//add coef to free_coefs
		else free_coefs[i] = free_coefs[i] + coef;
	      }
	      /*for(subtree_row_it sr_it = sr.begin(); sr_it != sr.end(); ++sr_it) {
		std::cout << "ASSUM VAR : " << *sr_it->first << " ";
		for(std::vector<double>::iterator i = (sr_it->second).begin();
		    i != (sr_it->second).end(); ++i) {
		  std::cout << *i << " ";
		}
		std::cout << std::endl;
		}*/
	    }
	    else { //no id::plus
	      //identify paires coef*expression
	      pre_it expression = assum_child;
	      contin_t coef =
		splitCoefExpression(ctr, expression) * assum_coef;
	      //add or uptdate (expression, rows) in the map
	      if(ctr.is_valid(expression)) {
		subtree_row_it sr_it = sr.find(expression);
		if(sr_it==sr.end()) {
		  std::vector<contin_t> v(assumptions.size()+2, 0.0);
		  v[i] = coef;
		  std::pair<pre_it, std::vector<contin_t> > p(expression, v);
		  sr.insert(p);
		}
		else { //update coef
		  std::vector<contin_t>& v_sr = sr_it->second;
		  v_sr[i] = v_sr[i] + coef;
		}
	      }
	      //add coef to free_coefs
	      else free_coefs[i] = free_coefs[i] + coef;
	    }  
	  }
	  //generate matrix to be solved
	  std::vector< std::vector<contin_t> > mat;
	  for(subtree_row_it sci = sr.begin(); sci != sr.end(); ++sci)
	    mat.push_back(sci->second);
	  mat.push_back(free_coefs);
	  //print matrix before
	  /*std::cout << "MATRIX BEFORE :" << std::endl;
	  for(unsigned int i = 0; i < mat.size(); i++) {
	    for(unsigned int j = 0; j < mat[i].size(); j++)
	      std::cout << mat[i][j] << " ";
	    std::cout << std::endl;
	    }*/
	  //perform Gaussian elimination, if fail return
	  bool resGJ = gaussianElimination2(mat);
	  //print matrix after
	  /*std::cout << "MATRIX AFTER :" << std::endl;
	  for(unsigned int i = 0; i < mat.size(); i++) {
	    for(unsigned int j = 0; j < mat[i].size(); j++)
	      std::cout << mat[i][j] << " ";
	    std::cout << std::endl;
	    }*/
	  //non linear decomposition
	  if(!resGJ) {
	    return;
	  }
	  //check if the linear decomposition is positive or negative
	  std::vector<contin_t> solution;
	  bool BSres = backSubstitution(mat, solution);
	  //print vector solution
	  /*std::cout << "SOLUTION :" << std::endl;
	  for(unsigned int i = 0; i < solution.size(); i++)
	    std::cout << solution[i] << " ";
	  std::cout << std::endl;
	  //print vector isStrictPositives
	  std::cout << "ISSTRICTPOSITIVES :" << std::endl;
	  for(unsigned int i = 0; i < isStrictPositives.size(); i++)
	    std::cout << (bool)isStrictPositives[i] << " ";
	    std::cout << std::endl;*/
	  if(!BSres) {
	    return;
	  }
	  bool existsStrictPositive = false;
	  bool existsStrictFromStrictPositive = false;
	  bool allPositiveOrNull = true;
	  for(unsigned int i = 0; i < solution.size(); i++) {
	    if(isStrictPositives[i]) {
	      if(solution[i] > 0) {
		existsStrictPositive = true;
		existsStrictFromStrictPositive = true;
	      }
	      else if(solution[i] < 0) {
		allPositiveOrNull = false;
	      }
	    }
	    else {
	      if(solution[i] < 0)
		existsStrictPositive = true;
	      else if(solution[i] > 0)
		allPositiveOrNull = false;
	    }
	  }
	  /*std::cout << "existsStrictPositive : " << existsStrictPositive << std::endl;
	  std::cout << "existsStrictFromStrictPositive : " << existsStrictFromStrictPositive << std::endl;
	  std::cout << "allPositiveOrNull : " << allPositiveOrNull << std::endl;*/
	  //positive reduction
	  if(allPositiveOrNull && existsStrictFromStrictPositive) {
	    //std::cout << "POSITIVE REDUCTION" << std::endl;
	    *it = (termStrictPositive? id::logical_true : id::logical_false);
	    tr.erase_children(it);
	  }
	  //reduction to false
	  else if(!existsStrictPositive) {
	    //std::cout << "NEGATIVE REDUCTION" << std::endl;
	    *it = (termStrictPositive? id::logical_false : id::logical_true);
	    tr.erase_children(it);	      
	  }
	}
      }
    }
    //std::cout << "LINEAR AFTER : " << tr << " IT : " << *it << std::endl;
  }

} //~namespace reduct

